Building Self- Consistent Triaxial Galaxy Models Using 

Schwarzschild's Method 
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ABSTRACT 

We use Schwarzschild's orbit superposition method to build self-consistent models of ellipti- 
cal galaxies with scale-free potentials. Our exhaustive study of all physical shapes and central 
densities for galaxies with scale-free potentials establishes a relationship between the presence of 
chaos and self-consistency. The extent and the onset of chaos is studied by varying parameters 
of the model which are believed to be its inducers, such as the steepness of the central density 
cusp, flatness and triaxiality. We show that gravitational scattering of the central density cusp 
plays a dominant role in restricting the shapes of elliptical galaxies. 



Subject headings: celestial mechanics 
ies: elliptical and lenticular, CD 

Introduction 



■ stellar dynamics — galaxies: kinematics and dynamics — galax- 



Recent observational data has shown that most, 
if not all, elliptical galaxies feature power-law den- 
sity distributions which rise into the smallest ob- 
servable radii, thus causing them to have central 
density cusps (Crane et al. 1993; Moller, Stiavelli 
& Zeilinger 1995; Lauer et al. 1995; Gebhardt 
et al. 1996). Inner parts of most of these galax- 
ies have two power-law regimes, with the core ra- 
dius being the boundary between the steep outer 
and shallow inner density profiles (Gebhardt et al. 
1996). Some, however, seem to obey one uniform 
power-law throughout, thus rendering them scale- 
free. It is also important to note that orbits in the 
outer parts of galaxies with two power-law regimes 
as well as the ones quite close to the center are es- 
sentially scale-free because of their distance from 
the core. These reasons, along with the fact that a 
great deal of simplification is introduced in study- 
ing this family of potentials, make scale-free mod- 
els of great interest. 

It has long been speculated that a presence of a 
black hole or a central density cusps may disrupt 
most box orbits and thus cause a loss of triaxi- 



ality in the innermost parts (Gerhard & Binney 
1985). Our goal is to study the dependence of 
the model's self-consistency on the strength of its 
central density cusp, as well as its flattening and 
triaxiality. In order to foster a thorough and com- 
prehensive investigation of this dependence, our 
models span the entire ranges of elongation and 
central density observed in nature, while varying 
triaxiality from one physical extreme to the other. 
Photometric observations of elliptical galaxies re- 
veal that they range from spherical EO to quite 
elongated E6, whose ratio of lengths of short to 
long axis is c/a — 0.4. High resolution HST ob- 
servations, which peer into the very centers of el- 
liptical galaxies, reveal that the galactic density 
profiles rise as power-laws p ~ r~ 7 as the center is 
approached, with the strength of the central cusp 
being in the range 7 = [0.25, 2] (Lauer et al. 1995). 

Schwarzschild's method has been a traditional 
tool for probing the dynamics of elliptical galax- 
ies. It was applied to a number of axisymmetric 
and disky (Cretton et al. 1999; Verolme & de 
Zeeuw 2002; Jalali & de Zeeuw 2002), as well as 
triaxial potentials, both integrable (Statler 1987; 
Siopis 1999) and non-integrable (Schwarzschild 
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1979, 1993; Merritt & Fridman 1996; Siopis 1999). 
Earlier applications of Schwarzschild's method to 
triaxial systems either studied models of one or 
two particular density profiles and varying galac- 
tic shapes (Statler 1987; Schwarzschild 1993; Mer- 
ritt 1997; Poon & Merritt 2001), or have investi- 
gated a range of density profiles for fixed galac- 
tic shapes (Merritt & Fridman 1996; Siopis 1999). 
Our study provides the most exhaustive coverage 
of both density profiles and galactic shapes, span- 
ning virtually the entire range of physically plau- 
sible models for galaxies with scale-free potentials 
and singular central densities. This enables us to 
better understand the dependence of models' self- 
consistency on these properties. By investigating 
orbits in scale-free potentials, we single out the ef- 
fects of power-law central density singularities on 
their stability. Our models present a useful tool for 
studying the stability and self-consistency of tri- 
axial elliptical galaxies with central density cusps. 
In §2, we outline the extent of our study and re- 
view major properties of scale-free potentials. §3 
covers major issues in the numerical implementa- 
tion of Schwarzschild's method. Major findings of 
the study are reported in §4 and their significance 
and impact discussed in §5. 

2. Mass Models 

We investigate scale-free models for which the 
isodensity surfaces are similar concentric ellip- 
soids. Their densities are given by the power-law 
formula 

(1) 
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We take a, b and c to be the long, intermediate, 
and short axes respectively. To simplify the com- 
putations of forces, we represent the density by 
the double expansion 
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in the usual spherical coordinates of radial dis- 
tance r, azimuthal angle 0, and polar angle 9 (Bin- 
ney & Tremaine 1987). Here Pf™ is an associated 
Legendre function. The special form of this ex- 
pansion is due to the eight-fold symmetry of the 



density which is even in x, y, and z. The potential 
is found from Poisson's equation as 
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for 7 = 2. Except for the special case of 7 = 2 and 
Tji — 71 — 0, the coefficients of the two expansions 
are related by the formula 
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[(2- 7 )(3-7)-2n(2n+l)]. (6) 



The coefficients of the expansion for the density 
are given by 
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where S m o is a Kronecker delta which is unity for 
m = but zero otherwise. Even with M max = 
4, the truncated Fourier-Legendre expansions (4) 
and (5) for the potential are accurate to within 
10~ 4 , while expansion (3) for the density is ac- 
curate to within 1% even for the most flattened 
models. 

3. Numerical Implementation 

Our computations made use of the IBM RS /6000 
SP supercomputer at the Florida State University 
School of Computational Science and Information 
Technology. It features 42 4-way 375 MHz Power 3 
nodes. Parallel programming and aggressive com- 
piler optimization resulted in twenty- to thirty-fold 
increases in speed from the more conventional se- 
rial runs on Ultra Sparc UNIX workstations. 

Orbit integrations used the DOP853 integra- 
tion routine, the explicit imbedded (7,8) pair of 
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Dormand and Prince (Hairer, Norsett & Wanner 
1993). It proved to be quite efficient in generat- 
ing large number of orbital points needed for the 
frequency analysis. The errors in the energy con- 
servation as well as the numerical integration have 
been on the order of 10~ 8 even after hundreds of 
dynamical times. 

3.1. Sampling the Phase Space 

In order to assemble a library of orbits which 
is representative of all orbits in a given poten- 
tial and without missing any major orbital fam- 
ilies, the phase space must be systematically sam- 
pled. Schwarzschild (1993) proposed a two-fold 
2-D start-space; the stationary start-space con- 
tains initial conditions starting from equipotential 
surfaces with zero velocities, while the principal- 
plane start-space features radially stratified initial 
conditions which pierce one of the three princi- 
pal planes with the velocity vector normal to the 
plane. These spaces are designed to pick up dif- 
ferent kinds of orbits arising in triaxial potentials; 
stationary start-space picks up orbits which have 
zero-velocity turning points, such as boxes and 
resonant boxlets, while the principal-plane start- 
space selects mostly tube orbits. It is probable, 
yet not certain, that the combination of these two 
start-spaces covers all the phase space of triaxial 
potentials (Schwarzschild 1993). In any case, it is 
necessary to understand which orbital structures 
are represented in such spaces, as well as deter- 
mine which ones are likely to be left out (cf. Ap- 
pendix A). 

We integrate 2304 orbits from the stationary 
start-space, by placing 12 equally spaced initial 
conditions in each of the 192 cells of the refer- 
ence sphere (cf. §3.4.1). A total of 1000 orbits 
with the same energy are sampled on each princi- 
pal plane, spanning an elliptical annulus with an 
outer radius equal to the radius of the equipoten- 
tial surface. The inner boundary of the annulus in 
each of the planes is chosen to be an ellipse whose 
axes are radii of the periodic 1:1 thin tube orbits 
perpendicular to that plane. This ensures that 
the innermost initial conditions from the principal 
plane start space correspond to thin tubes. In all, 
a total of 5304 initial conditions is integrated for 
each model (Figure 1). Each orbit is integrated 
for 200 dynamical times, which we define to be 
the maximum of the periods of thin tubes in three 



principal planes. 

a) z 




Fig. 1. — Initial condition spaces: a) stationary 
start-space b) principal-plane start-space. 



3.2. Orbital Structure 

We classify orbits based on the resonances be- 
tween leading frequencies in the Cartesian coor- 
dinates (LFCCs), and follow the paradigm given 
in the Table 1 of (Schwarzschild 1993). The most 
dominant families are boxes, short axis tubes (S- 
tubes), inner long axis tubes (I-tubes), outer long 
axis tubes (O-tubes) and resonant boxlets. The 
LFCCs of the boxes are all distinct; tube orbits, 
however, have the same LFCCs in coordinates nor- 
mal to their axis of rotation. For the resonant 
orbits, the ratio of LFCCs in two or more coordi- 
nates is a ratio of integers. Low order resonances, 
the ones for which the ratio of LFCCs is a ratio 
of small integers, occupy a larger portion of the 
phase space than the high-order resonant families 
(cf. Appendix A). This classification enables us 
to utilize the ratio of LFCCs space in analyzing 
populations of orbits. 

Two initial condition spaces provide two dis- 
tinct orbital populations. Stationary start-space 
yields mainly resonant orbits (boxlets) or unstable 
boxes, while the principal-plane start-space pro- 
duces mostly tubes. This duality is evident in 
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the graphs of the ratios of LFCCs. Note that the 
fundamental frequencies do not necessarily corre- 
spond to LFCCs; the two are equivalent for box or- 
bits and other resonant orbits, but the tube orbits 
have three distinct fundamental frequencies, while 
having the same LFCCs normal to the axis of rev- 
olution. This loss of information does not hurt our 
purposes; using LFCCs to classify orbits, measure 
chaotic diffusion and demonstrate duality of the 
two initial condition spaces is equivalent and just 
as effective as using fundamental frequencies. Fig- 
ure 2 shows the locations of orbits for each of the 
two start-spaces in the ratio of frequencies space. 
In the ratio of frequencies space, tubes are con- 
centrated around lines v x jv z — v y jv z (S-tubes) 
and v y jv z = 1 (I- and O-tubcs), which is indeed 
where most of the orbits from the principal-plane 
start-space are concentrated. Resonant orbits and 
boxes occupy a region between these lines, where 
a number of resonant lines intersect. Some overlap 
between the two spaces is to be expected, since nei- 
ther can exclusively pick up only one population of 
orbits. It is particularly interesting to look at the 
transition from prolate to oblate models. For pro- 
late models, all orbits are I- or O-tubes, located 
along the line v y jv z — 1, but as the triaxiality 
decreases and the models become more oblate, or- 
bits in both spaces seem to shift steadily toward 
the S-tubes (line v x jv z = v y /v z ) until they com- 
pletely align with it for oblate models. Our earlier 
studies of 2D oblate and prolate scale-free mod- 
els revealed very little chaos (Hunter et al. 1998), 
which is consistent with Figure 2. The bottom row 
shows the fraction of regular orbits as a function 
of the chaotic diffusion parameter ui, as defined in 
(9). 

3.3. Detecting Chaotic Orbits 

Chaotic orbits are not time-independent - their 
orbital properties, such as eccentricity, radial 
excursion and orbital densities, evolve in time. 
In three degrees of freedom systems, all chaotic 
portions of the phase space are interconnected 
(Arnold's web), so that every chaotic orbit will 
eventually visit each chaotic section. Thus, we 
may choose to look at each of the chaotic or- 
bits as different parts of one chaotic 'super-orbit'. 
This super-orbit gives the averaged density of the 
stochastic portions of the phase space, and there- 
fore is time- independent. If orbit integration is 
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Fig. 2. — Ratios of leading frequencies in Carte- 
sian coordinates for the stationary start-space 
(first column) and principal-plane start-space (sec- 
ond column) for (top to bottom) T=l (prolate), 
T=0.7 T=0.3 and T=0.0 (oblate). The third col- 
umn shows the fraction of orbits in each of the 
four major orbital families: B - boxes, S - short 
axis tubes, L - long axis tubes, R - resonant or- 
bits. The fourth column shows the fraction of or- 
bits with the diffusion parameter greater than u. 
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carried over sufficiently long time intervals, then 
the time-averaged orbital density of each chaotic 
orbit would approximate the density of the chaotic 
super-orbit. However, the integration interval re- 
quired for achieving a good time-averaged approx- 
imation may be quite long for some weakly chaotic 
orbits residing near the boundary between chaotic 
and regular portions of the phase space, thus ren- 
dering this approach numerically impractical. It 
is more efficient to take a simple arithmetic aver- 
age of orbital densities of all chaotic orbits of the 
model to get a more accurate representation of 
the chaotic super-orbit (Merritt & Fridman 1996). 
We are looking for true equilibrium solutions, and 
therefore only include time independent building 
blocks into Schwarzschild's method: regular orbits 
and one chaotic super-orbit. This necessitates 
careful distinction between chaotic and regular 
orbits. 

Both regular orbits and the chaotic super-orbit 
can be viewed as ergodic - they sample all of the 
allowed phase space. A major difference between 
the two is that the phase-space volume occupied 
by the chaotic orbits is much larger, thus making it 
impossible to represent it accurately with just one 
orbital 'snapshot' obtained after a finite integra- 
tion of a single orbit. The size of the phase-space 
volume of the stochastic region enables vast num- 
bers of these different orbital snapshots of a same 
super-orbit to exist. Therefore, in order to have 
a good approximation to the contribution of the 
stochastic orbits to the model's density, we com- 
bine the information conveyed in the individual 
snapshots into a single chaotic super-orbit. 

We use Hunter's method (Hunter 2002) to 
compute the diffusion in fundamental frequencies 
(FFs) of an orbit integrated over two consecutive 
time intervals. The method extracts frequencies 
from an orbit given in the form of a time series. 
We have discovered that Hunter's method is more 
accurate and more efficient than Laskar's numer- 
ical analysis of fundamental frequencies method 
(Laskar et al. 1992; Laskar 1993) by a few orders 
of magnitude. This is because it relies entirely on 
discrete Fourier transforms. All of the extracted 
frequencies will be linear combinations of three 
FFs. These FFs will remain constant if the mo- 
tion is regular. When the motion is chaotic, these 
FFs will change. The relative rate of change in 
FFs of an orbit over the integration interval is 



the measure of its averaged stochasticity over that 
time span. We compute the norm of the diffusion 
in the leading frequencies in Cartesian coordi- 
nates, instead of the norm of the diffusion in the 
fundamental frequencies. The two are equivalent, 
since frequencies in Cartesian coordinates are lin- 
ear combinations of the fundamental frequencies. 

We seek to compute the average stochasticity of 
an orbit over a certain finite time interval. The av- 
erage rate of diffusion in FFs is computed for each 
orbit over an integration interval t = [0, 2T] , where 
2T corresponds to nd = 200 dynamical times. We 
use these rates of diffusion to compute how many 
crossings it would, at that rate, take for a partic- 
ular orbit to 'lose memory' of its initial FFs, i.e. 
become relaxed. More precisely, if we measure the 
FFs over two consecutive time intervals t\ = [0, T] 
and t 2 — [T, 2T] to be tti and fi 2 , the average 
diffusion rate of FFs over the time interval [0, 2T] 
is 

m - ISi^kl. (8) 

We define an orbit to be relaxed when the change 
in FFs in on the order of the FFs themselves. At 
the rate of diffusion of SCI, this corresponds to the 
number 

of crossings (or dynamical times). Measuring dif- 
fusion in FFs is equivalent to computing short time 
Lyapunov exponents, since both measure an or- 
bit's average orbital stochasticity over a finite time 
interval. 

The scale-free potentials do not have an intrin- 
sic time/length scale associated with them. This 
requires our criterion for determining an orbit's 
stochasticity to be tied in with the number of 
crossings of an orbit instead of some absolute time 
interval, such as the Hubble time. We define u) to 
be the number of crossings it takes for an orbit 
to become fully relaxed. We adopt a convention 
that an orbit is chaotic if the full relaxation is ex- 



pected to set in within 10 4 crossings (\og 1 



4). 



Considering that the number of crossings a typi- 
cal orbit of an elliptical galaxy has gone through is 
between 20-50 for the ones in the outer parts and 
a few hundred for the ones near the center, this 
seems like a reasonable cut-off point. We discuss 
the effects of varying this stochasticity criterion in 
§4.1. 
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As the chaotic threshold is increased, more or- 
bits are labeled as chaotic and included into a 
super-orbit. In theory, averaging over all chaotic 
orbits will approximate the phase-space density of 
the interconnected chaotic sea, with its accuracy 
increasing as with the increase in the number of 
stochastic orbits and the length of the integration 
time. In practice, however, the integration time 
is finite, and the weakly chaotic orbits may only 
sample a limited region of the phase space, which 
may lead to a super-orbit which is not entirely 
time- independent. Possible time-dependence of 
the super-orbit is of little practical significance, 
since we find that including the super-orbit, as 
opposed to not including it, will never change 
model's self-consistency. This is the consequence 
of the fact that the super-orbit is an average of 
an ensemble of many nearly-ergodic orbital den- 
sities, which make its orbital density round and, 
as such, not crucial for reproducing the desired 
galaxy shapes. 

3.4. Constructing Self-consistent Scale- 
Free Models 

3.4-1- Coarse-graining of the Configuration Space 

The advantage of using a scale-free potential 
with triaxial symmetry is that one needs to con- 
sider only an octant of a 2-D reference sphere 
(r = R re{ ) (Richstonc 1980; Schwarzschild 1993). 
Each orbit produces a template orbital density 
which represents an ensemble of geometrically 
similar orbits. For any orbit through a point 
at radial distance r, there will be a geometri- 
cally similar orbit which passes through the refer- 
ence sphere whose length scale differs by a factor 
R Te f/r. With density proportional to r~ 7 , it fol- 
lows that we should give each orbital point on a 
computed orbit the weight (r/i? rc f) 7-3 (Richstone 
1980; Schwarzschild 1993). 

Following Schwarzschild (1993), we consider 
only a positive octant of the reference sphere, 
which we divide into 3 regions with planes x = y, 
x = z, and y — z, so that Region 1 is the one 
with x > y, x > z, Region 2 with y > x, y > z 
and Region 3 with z > x, z > y. Each region 
is then subdivided into 64 parts of equal surface 
area as follows. We first divide Region 1 into two 
equal halves with a vertical line v = y/x = Vi, 
and then split each half in two with two horizon- 
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Table 1 : Limits of integration for each of the cells 
on the reference sphere. 

tal lines h = z/x — li and h — ri. This step 
yields four parts of equal surface area, and can be 
recursively applied to each partition; in n such re- 
cursive steps, the region can be divided into 4" 
cells of equal area. We repeat this process un- 
til the region is split in 64 equal parts. Similar 
partitioning is done for the other two regions by 
rotating coordinates (x,y,z) — > (y,x,z) for Re- 
gion 2, and (x, y, z) — > (z, x, y) for Region 3. We 
use a Maple program to carry out the relatively 
straightforward integral computations involved in 
finding the positions of the vertical and horizon- 
tal dividing lines, i.e. the values of Vi, li and r^. 
Figure 3 shows the entire equipartitioned octant. 
Schwarzschild (1993) uses cells that are different 
in area to within few percent. Even though the 
cquipartitioning of the reference sphere is not one 
of the key issues in the implementation of the 
Schwarzschild's method, it is, nonetheless, note- 
worthy that the algorithm outlined above achieves 
maximal accuracy with the same amount of work. 

3.4-2. Computing Model Mass in Cells 

We compute surface mass on the reference 
sphere by integrating the density (3) over r = 
R rc f = 1 to get 
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Fig. 3. — Equipartitioning of the reference sphere. 



for the surface mass in the i-th cell. Each cell on 
the reference sphere is delimited by two constant 
values of each horizontal and vertical coordinates 
(h, v), which, in turn, provide limits of integration 
for <f> and 8 (Table 1). The inner integrand is a 
polynomial of degree 2n in cos#, for which the 
Gaussian quadrature with 2n points is exact. 

3.4-3- Computing Orbital Mass in Cells 

We need to be able to compute accurately the 
amount of time orbits spend in each of the cells 
of the reference sphere. When integrated for a 
fixed time interval, an orbit will generally not go 
through an integer number of crossings. In par- 
ticular, if an orbit has n complete crossings in t, 
but t/n is not an integer, some parts of the orbit's 
path will be sampled n + 1 times, while others 
only n times. This will introduce a certain inaccu- 
racy, which we minimize by increasing t in order to 
make n large. Orbits are given by a dense output 
of orbital points, which we further improve upon 
by linearly interpolating between them. We use 
an algorithm proposed by Siopis (1999); as long 
as consecutive orbital points are inside the same 
cell, advance orbital output in steps of to. When 



the next orbital point is outside the current cell, 
the time is advanced in steps that are a factor of 
10 and 100 smaller as the boundary between cells 
is approached. This algorithm ensures that the ac- 
curacy in determining the orbit's contribution to 
the mass of each cell is computed to within 0.01 to- 

3.4-4- Formulating the Optimization Problem 

Schwarzschild's method is formulated as an op- 
timization problem: 

Minimize : / (wi), 
Subject to : w i Pij = Ph (H) 

»=1 

Wi > 0, 

where j = 1, 2, N c , i = 1, 2, N Q , f(wi) is the 
cost function, is the contribution of the or- 
bital density of the ith template to jth cell, pj 
is the model's density in the jth cell and Wi is 
the orbital weight of the ith orbit. This becomes 
a linear programming problem (LPP) when the 
cost function / is a linear function of the weights; 
for example, to minimize weights of orbits labeled 
from to to n, the cost function would simply be 

n 

f(wi) — w i- The solutions of the LPP are 

i—m 

often quite noisy, with a significant number of or- 
bits unpopulated (Siopis 1999). It is often cus- 
tomary to impose additional constraints in order 
to 'smoothen' out the solutions, such as minimiz- 
ing the sum of squares of orbital weights, which 
makes this a quadratic programming problem, or 
minimizing the least squares. A thorough discus- 
sion of effects of this 'smoothing', as well as prob- 
lems involving the formulation of the optimization 
problem, is given in Siopis (1999). Since our study 
focuses on establishing existence of self-consistent 
solutions, we solve only the LPP. We use the 
BPMPD (Meszaros 1996) routine, a primal-dual 
interior point algorithm for LPP. It features pow- 
erful pre-solve techniques, efficient sparcity han- 
dling, and flexible linear algebra. Each LPP of 
size N Q x 192, where N Q < 5304, takes only a few 
seconds of CPU time on the RS/6000 SP. 

In order to investigate how non-self-consistency 
sets in, it is advantageous to know which cells of 
the reference sphere are infeasible. We do this by 
reformulating the LPP above into another LPP 
which uses slack variables Xi and in to make the 
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problem always feasible (Siopis 1999): 

N c 

Minimize : / (wi) +p^^(Aj + 
3=1 

Subject to : \j — + ^ Wi pij = pj, (12) 

where j = 1,2, ...,N C , i — 1,2, ...,N a , and p, the 
penalty parameter, is some positive constant. A 
self-consistent problem will have zero Xj and fij, 
while the deviation from exact self-consistency in 
non-self-consistent models can be traced by locat- 
ing the non- vanishing values of the slack variables. 

4. Results 

We study the set 7 = {0.5,1,1.5,2} of values 
of the central density cusp. For each 7 we vary 
the elongation from c/a = 0.4 to 1 in increments 
of 0.1. At each elongation, the triaxiality is var- 
ied through the triaxiality parameter T = (a 2 — 
b 2 )/(a 2 - c 2 ), from oblate (T = 0) to prolate (T = 
1) in the interval T = {0,0.1,0.3,0.5,0.7,0.9,1}. 
Figure 4 plots the 42 different shapes we investi- 
gated in axis-ratio space. The hypotenuse of the 
triangle (line b — c) represents prolate spheroids 
(T = 1), the vertical side (b = a line) represents 
oblate spheroids (T = 0), the vertex (1,1) at which 
those two sides meet represents spherical models 
(Figure 4). 

Only regular orbits and a single chaotic super- 
orbit (as an approximation to a time-independent 
phase space density of the chaotic portion of the 
phase space) are used in building self-consistent 
models using Schwarzschild's method. Figure 5 
shows the fraction of orbits from both station- 
ary and principal-plane start-spaces, as a func- 
tion of the degree of stochasticity ui for models 
with weak and strong central cusps, 7 = 0.5 and 
7 = 2 respectively. We note that the curves rep- 
resenting each of the start-spaces are relatively 
close together. As the strength of the cusp in- 
creases, the two lines move in opposite directions 
for non-axisymmetric cases (T ^ 0,1); the one 
representing the orbits from the stationary start- 
space drops off sooner and more abruptly, while 
the line corresponding to the principal-plane start- 
space becomes rounder, with its drop-off delayed. 
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Fig. 4. — Axis ratio space: the stars represent 

models investigated. The lines connect points 

of equal triaxiality, from right to left, T = 
0,0.1,0.3,0.5,0.7,0.9,1. 



This clearly implies that increasing the strength 
of the central cusp has a destabilizing effect on or- 
bits from the stationary start-space, which, as we 
established earlier, are mainly boxes and boxlets 
which pass near to the center. Increasing central 
mass concentration makes gravitational scattering 
of orbits more efficient, thus causing orbits to be- 
come more chaotic. Also, the tube orbits, which 
dominate the principal-plane start-space, are sta- 
bilized by the strengthening of the central cusp. 
They stay significantly away from the center, and 
therefore 'see' the central cusp as a point mass. 
The strengthening of the central cusp has the same 
effect on them as increasing the central point mass 
- it moves them toward a more Keplerian, inte- 
grable regime, causing them to be less chaotic. 
This behavior is more pronounced for the rounder 
models. In the absence of a central density cusp, 
we expect rounder, more spherical, models to be 
closer to integrability and thus to feature more reg- 
ular orbits. This in turn means that the chaotic 
effects of introducing a central density cusp are 
better isolated and therefore more prominently 
displayed in rounder models. In the case of ax- 
isymmetric models, most orbits are regular tubes 



8 



whose stability is only reinforced by the central 
mass concentration. 
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Fig. 5. — Fraction of orbits as a function of 
the degree of stochasticity for models: a) 7 = 
0.5, b) 7 = 2. Dashed lines represents orbits from 
the stationary start-spaces, and the full lines orbits 
from the principal-plane start-spaces. log 10 u> = 4 
is the chaotic threshold. Each column represents 
models with triaxiality T = 1.0,0.7,0.3,0.0 (from 
left to right) and each row models with long-to- 
short axis ratios c/a — 0.9,0.7,0.4 (top to bot- 
tom). 

In Figure 6, we analyze the orbital structure of 
models by graphing the fraction of orbits in each 
of the four major orbital families: boxes, short 
axis tubes, long axis tubes and low-order reso- 
nant orbits. One drawback of this coarse break- 
down is that higher order resonant orbits will be 
identified as boxes. As expected, bounding ax- 
isymmetric models will feature mostly short axis 
tubes (oblate, T = 0) and long axis tubes (prolate, 
T = 1). As the central density cusp is increased, 



a roughly even distribution between the boxes, 
short axis tubes and long axis tubes, present in 
weak-cusped models only for more spherical cases, 
expands to include more flattened models with 
steeper central density cusps. If we recall that box- 
dominated stationary start-space composes only 
about 40% of the total number of orbits inte- 
grated, it becomes evident that many initial condi- 
tions from the principal-plane start-spaces of flat- 
ter, weak-cusped models also produce box orbits. 
The work of Statler (Statler 1987) on building 
self-consistent models of 'perfect triaxial galaxies' 
which feature smooth cores has evidenced similar 
behavior to our weak-cusped case; as one moves 
down the middle of the figure but away from the 
oblate and prolate limits, the box orbits become 
progressively more dominant. The reason for this 
is that, as the mass model becomes elongated in 
one direction, the tube orbits elongate in the di- 
rection normal to it, thus opposing the shape of 
the mass distribution. This leaves the box orbits 
as the only ones which can support the triaxial 
shapes for most models (Statler 1987). 

The solutions of the LPP formulation of the 
Schwarzschild's method are summarized in Figure 
7. For weak cusps, chaotic scattering is less effi- 
cient and many regular orbits suitable for repro- 
ducing broad ranges of galactic shapes still escape 
the destabilizing effects of the center. This re- 
sults in finding self-consistent solutions for all but 
the most flattened shapes of galaxies with weak 
central density cusps (Figure 7. a). The non-self- 
consistency first sets in for flat and nearly pro- 
late models. As the strength of the central cusp 
is increased, the region of non-self-consistency ex- 
pands along the prolate boundary (without actu- 
ally including it) and up the middle of the axis 
ratio space to include more triaxial models. Even- 
tually, in the case of the strong central density 
cusp, i.e. the logarithmic scale-free potential, the 
region of non-self-consistency dominates the en- 
tire axis ratio space. It is important to note that 
the self-consistent region around the oblate edge 
of the axis ratio space remains fairly large, even 
when a strong central density cusp is present, 
while the one around the prolate models remains 
very thin. This is consistent with the previous 
findings in which self-consistency prevails only for 
nearly oblate and prolate models; Figure 1 of 
(Merritt 1997) shows that for a Jaffe's density 
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power law, which behaves as p ~ r~ 2 near the 
center, i.e. has a strong central density cusp of 
7 = 2, the self-consistent models are confined to 
a narrow strip near the prolate boundary and a 
much thicker region around the oblate boundary. 
Hollcy-Bockclmann ct al. (2001) used TV-body 
simulations to construct self-consistent models for 
a Hernquist 7 = 1 profile with a range of triaxial- 
ity and modest flattening, which is consistent with 
our Figure 7.b. 
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Fig. 7. — Self-consistent (•) and non-self- 
consistent (o) models for each of the four models: 
a) 7 = 0.5, b)7=l, c) 7 = 1.5, d) 7 = 2. 



Fig. 6. — Fraction of orbits in each of the five 
orbital families: B - boxes, S - short axis tubes, 
L - long axis tubes, R - second and third or- 
der resonant orbits, a) 7 = 0.5, b) 7 = 2. 
Each column represents models with triaxiality 
T = 1.0,0.7,0.3,0.0 (from left to right) and each 
row models with long-to-short axis ratios c/a = 
0.9,0.7,0.4 (top to bottom). 



In order to understand the transition from self- 
consistency to non-self-consistency we look at the 
infcasiblc cells of the reference sphere, the ones for 
which constraints (12) are not satisfied. Figure 8 
shows the non-feasible cells for each of the models 
(denoted by dots). As the strength of the cen- 
tral density cusp increases, non-self-consistency 
is introduced through infeasibility of cells mainly 
around the short axis (z-axis) and the y-z plane. 
These infeasible regions correspond to chaotic re- 
gions of the stationary start-space. Orbits which 
would reproduce the model's density in those in- 
feasible regions are chaotic and as such averaged 
into a round super-orbit. We solidify this argu- 
ment by integrating additional dense population of 
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orbits in these infeasible regions, which results in 
most orbits being chaotic, while the regular ones 
were still not sufficient to reproduce the desired 
density. This strongly suggests that the chaotic 
nature of orbits needed to reproduce mass den- 
sity in those infeasible cells is responsible for the 
non-self-consistency. However, before any mean- 
ingful conclusions about the intrinsic dynamics 
can be made, one must clearly understand the im- 
plications and limitations of the Schwarzschild's 
method and its dependence on parameters such 
as the number of orbits integrated, the number 
of cells of the reference sphere, and the stochastic 
threshold. 

a) T=0.9 T=0.7 T=0.5 T=0.3 T=0.1 

c /a=o, 8 ' , : a a r\ / 3 
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Fig. 8. — Non-feasible cells (marked by dots) for 
models with a) 7 = 1 and b) 7 = 1.5 on a grid 
of 192 cells. 



4.1. Distinguishing Between Numerical 
Artifacts and Intrinsic Dynamics of 
the System 

As with any numerical model, it is necessary 
to filter out all of the numerical effects associ- 
ated with Schwarzschild's method from the un- 
derlying dynamics of the system before any phys- 
ical conclusions are reached. There are two ma- 
jor assumptions and simplifications introduced by 
Schwarzschild's method. First, the configuration 
space is coarse-grained, divided into cells on which 
we require the distribution function (DF) to be 
self-consistent, i.e. to satisfy both the collision- 
less Boltzmann equation (CBE) and the Poisson 
equation simultaneously. Second, the system is 
in equilibrium, i.e. the phase space DF is time- 
independent. 

We probe the coarse-grainedness of the solution 
by refining the reference grid - we test whether 
the solutions which are self-consistent on the orig- 
inal grid of 192 cells remain self-consistent when 
the grid resolution is doubled. The LPP is solved 
for this finer grid of 384 cells for the 7 = 1.5 
model (Figure 9; compare to Figure 8.b). Only 
oblate and prolate models remain self-consistent, 
while other previously self-consistent models fail 
to satisfy constraints of the optimization prob- 
lem in several cells. Mathematically, the num- 
ber of constraints increases with the number of 
cells, which shrinks the solution space of the LPP 
and eventually renders the problem infeasible. By 
the same token, increasing the number of orbital 
templates increases the number of free variables 
in the LPP and thus expands the solution space. 
For non-self-consistent models with 192 cells, in- 
creasing the number of orbital templates will not 
yield self-consistent solutions, which shows that 
the solutions above reflect true dynamical prop- 
erties of the model, rather than a numerical arti- 
fact. On the other hand, models which are self- 
consistent on a grid of 192 cells will remain self- 
consistent if both the number of orbital templates 
and cells (free variables and constraints, respec- 
tively) are increased simultaneously. This strongly 
implies the existence of self-consistent solutions in 
the continuous limit. We can, therefore, study the 
effects of central density cusps and galaxy shapes 
on the true self-consistency (the continuous limit) 
by studying self-consistency on the discrete grid of 
192 cells. 
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Fig. 9. — Non-feasible cells (marked by dots) for 
7 = 1.5 on a grid of 384 cells. 

The integration time of 200 dynamical times 
was chosen so that the orbital density is well sam- 
pled. Such long integration times ensure that most 
of the 'sticky' orbits (Contopoulos 1971; Kan- 
drup, Pogorelov & Sideris 2000; Siopis & Kan- 
drup 2000), which behave as regular for many 
dynamical times, are detected. If the orbit inte- 
gration time was any shorter, and the stochastic- 
ity criterion not stringent enough to detect such 
weak stochasticity, many of these orbits would 
have been included in Schwarzschild's method as 
regular, time-independent building blocks. We 
find that if both regular and chaotic orbits are 
included into Schwarzschild's method as time- 
independent building blocks, self-consistent solu- 
tions are readily found for most models (Figure 
10, first column). These are non-equilibrium so- 
lutions, producing a time-dependent phase space 
distribution function. Schwarzschild (1993) built 
several flattened, strongly triaxial, self-consistent 
models which included chaotic orbits integrated 
over one Hubble time. Upon integrating orbits 
from his self-consistent models over three Hub- 
ble times instead, the change in orbital densities 
of the stochastic orbits caused the models to be- 
come more round. However, including the time- 
dependent building blocks in the search for time- 
independent DFs leads to time-dependent solu- 
tions, and violates the initial assumption that the 
system is in equilibrium. If a numerically ob- 



tained DF, a superposition of orbital properties of 
the library of integrated orbits, undergoes albeit 
secular changes, the resulting change in model's 
mass density gives rise to a different set of or- 
bits. Then the previously computed DF is not a 
result of a superposition of orbits arising from the 
changed mass density, and therefore will not be a 
true self-consistent solution, i.e. will not solve Poi- 
son's equation and CBE simultaneously. However, 
in practice it is useful to relax the strict equilib- 
rium requirement, in order to allow those slowly- 
evolving, quasi-equilibrium solutions to provide us 
with an insight as to what may happen with self- 
consistent models in time. Their steady evolution 
toward more rounder shapes, leads one to suspect 
that the triaxiality may be a transient feature of 
elliptical galaxies. 



chaotic orbits logio (0=3 logio CO = 4 logio CO = 5 
included , 




Fig. 10. — Infeasible regions of the axis ratio space 
(shaded) when chaotic orbits are treated as regular 
(first column) and for varying values of chaotic 
threshold ui. 



It is interesting to observe the dependence of 
self-consistent solutions on varying the chaotic cri- 
terion used to distinguish between regular and 
chaotic orbits. As the chaotic criterion becomes 
less restrictive, i.e. as the stochastic threshold 
u) decreases, a larger number of orbital templates 
are used in the search for self-consistency. Fig- 
ure 10 demonstrates that as to is decreased, it is 
possible for a previously non-self-consistent model 
to attain self-consistency. That is indeed what 
happens for most non-self-consistent models under 
consideration. If all orbits are treated as regular 
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and included in the Schwarzschild's method, self- 
consistency is attained for all but the flattest mod- 
els with strong cusps. On the other hand, if only 
regular orbits are included (Merritt 1997), or if 
they are combined with a chaotic super-orbit, the 
self-consistent regions are significantly limited as 
we show in this study. This is the consequence of 
the fact that in flatter and more cuspy potentials, 
the regular orbits do not have sufficient variety, or, 
more precisely, do not reproduce the density of the 
model along the axes. The orbits that would sup- 
port the prescribed density in these potentials and 
thus foster self-consistency arc rendered chaotic by 
the scattering effects of the central mass. This 
points to chaos as the main cause of non-self- 
consistency. We also observe that as the central 
density cusp becomes steeper, more centrophilic 
orbits become chaotic through gravitational scat- 
tering of the massive center, thus causing more 
models to become non-self-consistent through the 
shortage of regular time- independent orbits which 
can reproduce the mass density of flattened and 
triaxial models. It is evident from this that the 
gravitational scattering by the massive center is 
the main inducer of chaos. 

It has long been known that the solutions to 
self-consistent problems are highly non-unique 
(Hunter 1995; Statler 1987; Merritt 1997). Em- 
pirically we see this from the fact that we were 
able to find as many different DFs as we had cost 
functions for the LPP problem. Our results show 
that whenever self-consistent solutions are found, 
there exists a solution which includes only the reg- 
ular orbits. Such a solution is found by choosing 
a cost function for the LPP such as to minimize 
the chaotic super-orbit. This is, indeed, to be 
expected, since the orbital density of the chaotic 
super-orbit is round, and not particularly suitable 
for reproducing mass density for elongated and 
strongly triaxial models. 

5. Discussion and Conclusion 

We have studied the effects of central density 
concentrations on the existence of self-consistent 
solutions computed using Schwarzschild's orbit su- 
perposition method. The ranges of central den- 
sity cusps, elongation, and triaxiality, render this 
analysis comprehensive in its scope. The imple- 
mentation of the method is described and justi- 



fied in detail. Some of its main aspects provide 
a great deal of information about the dynamics 
of the system. The start-spaces, stationary and 
the principal-plane, provide us with a systematic 
manner in which to sample the phase space and 
gain an insight about the orbital structure of a 
given triaxial potential. The classification of or- 
bits based on the leading frequencies in Cartesian 
coordinates provides us with further information 
about the orbital structure. Detecting chaotic or- 
bits using a frequency extraction algorithm based 
on discrete Fourier transforms enables us to distin- 
guish between stochastic and regular orbits. This 
is essential when it comes to including them in 
Schwarzschild's method, since the two are treated 
differently. We also analyze and discuss in detail 
the assumptions, implications, and limitations, of 
Schwarzschild's orbital superposition method. 

The self-consistent solutions obtained through 
Schwarzschild's method with chaotic orbits aver- 
aged into a super-orbit are equilibrium solutions, 
since all of their constituents are time- independent 
building blocks. This time-independence restric- 
tion is proven to be the reason why a number 
of flattened triaxial models are found to be non- 
self-consistent; some orbits necessary to reproduce 
model's elongation and triaxiality are chaotic and 
thus averaged into a more round super-orbit, caus- 
ing the mass constraint not to be satisfied. On the 
other hand, if chaotic orbits are directly included 
in Schwarzschild's method, irrespective of the fact 
that their orbital properties evolve in time, such 
quasi-equilibrium self-consistent solutions are eas- 
ily found for most models. This clearly implies 
that there is no shortage of orbits which could 
temporarily reinforce triaxiality even for time in- 
tervals several times longer than the age of the 
Universe. Truly permanent existence depends on 
the regularity of orbits responsible for reinforcing 
the prescribed galactic density distribution. 

The scale-free property of the potentials inves- 
tigated here restricts us from knowing directly 
anything about the self-consistency of models at 
varying time/length scales. However, it does en- 
able us to isolate the effects of the central density 
cusp on the global dynamics of the system. Our 
findings strongly suggest that strengthening the 
central density cusp increases its ability to scat- 
ter centrophilic orbits efficiently and render them 
chaotic. This is in agreement with earlier stud- 
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ies of cuspy potentials (Merritt & Fridman 1996; 
Merritt 1997; Siopis 1999). We also find that this 
results in increased stability of the centrophobic 
tube orbits, since for them the strengthening the 
central density cusp has the effect similar to that 
of the increase of the central point mass on Ke- 
plcrian orbits. The breadth of our study, span- 
ning virtually all plausible galactic shapes, elon- 
gations, and central densities, shows us how and 
to what extent triaxiality and central density can 
cusps coexist. We find self-consistent solutions of 
weakly-cusped galaxies for almost the entire range 
of triaxial shapes, while the self-consistent region 
of the axis-ratio space for strong cusps is limited to 
nearly axisymmetric, mildly triaxial, regions near 
the prolate and oblate boundaries. 

Our study shows that the gravitational scatter- 
ing of the massive galactic center is more effective 
in rendering centrophilic orbits chaotic as the cen- 
tral density cusp becomes stronger. The scattering 
rids models of regular box-like orbits necessary to 
reproduce flattened triaxial shapes, thus rendering 
them non-self-consistent. This establishes gravita- 
tional scattering as the key factor in restricting the 
shapes of elliptical galaxies. 

The author is grateful to Christopher Hunter 
for stimulating discussions, comments, and sug- 
gestions. Suggestions by the anonymous referee 
contributed to the clarity and accuracy of the 
manuscript. The Florida State University School 
of Computational Science and Information Tech- 
nology has been generous in granting access to 
their supercomputer facilities, which greatly ex- 
pedited the computations presented here. This 
study has been supported by the NSF grant DMS- 
9704615. 
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A. Orbital Structure of 2-D Start Spaces 

The stationary start-space samples only the positive octant of the equipotcntial surface. This poses 
no restrictions on the orbits represented because of the triaxial nature of the potential. Initial conditions 
are sampled in some systematic fashion over this octant: one can either select the vertices of the equally 
spaced rectangular grid in each of three regions bounded by planes x = y, x = z and y = z (Schwarzschild 
1993; Mcrritt & Fridman 1996) or choose centers of equal area segments (Siopis 1999). There is no obvious 
advantage of one choice over the other. 

The principal plane start-space was initially thought of as only an x-z start-space (Schwarzschild 1993), 
but was later extended to the other two principal planes so that the symmetries of the system are better 
respected (Siopis 1999). In both cases, an attempt to minimize duplication of orbits is made by restricting 
the portion of the principal plane sampled to an elliptical annulus having the equipotential surface as the 
outer, and the minimum of the amplitudes of the 1:1 periodic thin tube orbits perpendicular to that plane 
as the inner boundary. The justification of this choice is as follows. Each tube orbit which crosses a given 
principal plane does so at two points, except in the case of thin tubes, for which the two merge (Statler 
1987). Therefore, Schwarzschild argues (Schwarzschild 1993), if one determines the line in the start-space 
on which thin tubes are located, and if the areas inside these thin tube lines are ignored, one is left with 
the part of the start-space in which each tube or box orbit is represented by one point (Schwarzschild 1993). 
This indeed will rid us of the unwanted non-uniqueness of the tube orbits, but it will also eliminate some 
other non-tube orbits entirely. We illustrate this on the same example as Schwarzschild used to argue his 
point: an x-z start-space for the 2-D oblate logarithmic potential with c/a = 0.3 and T = (Schwarzschild 
1993, Figure 2). 

For 2-D potentials, the x-z start-space represents both stationary and the principal-plane start-space. 
Each point in the space is a zero- velocity turning point of an orbit at some L z , as well as a point at which 
an orbit pierces the principal plane with velocity component normal to that plane (if one recalls that the 
magnitude of the z-component of the angular momentum is given by L z — \xy — yx\). It may be viewed 
as a compression of the continuum of surfaces of section: one for each curve of constant z-component of 
the angular momentum (L z ) (Schwarzschild 1993). Each point on that curve represents an orbit at that 
angular momentum. A point at which the 1:1 thin tube orbit touches this curve is the boundary which 
Schwarzschild proposes for the inner boundary of the ellipsoidal annulus. This leaves out all orbits which 
touch the L z = const, curve only at points to the left of this point, such as the 2:3 'fish' orbit, 2:5 and other 
resonant orbits (Figure 11, top panel). These orbits occupy a non-negligible portion of the phase space, as 
can be seen from the traditional Poincare surfaces of section (Figure 11, bottom panel). It is evident from 
the bottom panel of Figure 1 1 that the order of the resonance of the boxlet orbit is inversely proportional 
to the area of the phase space that it occupies. 

In order to better understand the structure of orbits on a L z = const, line, we compute the ratios of 
leading frequencies in x and z coordinates for a large number of orbits with initial conditions equally spaced 
along the curve. Figure 12 shows such ratios for orbits along curves as L z is increased. All thick tube orbits 
are represented by two points, one on each side of the thin tube. Clearly, the ratio of frequencies at such 
two points is the same, since they both represent the same orbit. Therefore, if the phase space consists 
only of tubes, the graph of the ratios of leading frequencies should be symmetric with respect to the thin 
tube location (as in Figure 12. f), albeit with x-scale somewhat altered. Furthermore, any deviations from 
such symmetry indicate the presence of resonant non-tube orbits. These resonant orbits are manifested in 
the graphs of ratios of leading frequencies by sudden jumps and flattenings away from the thin tube. The 
intersection of the plane in which the Poincare surface of section is taken with the thick resonant orbit is a 
set of invariant 'islands' surrounding the thin resonant orbit. Our earlier study (Terzic 1998) showed that 
all of these islands have the same rotation number (defined to be the fraction of a full revolution that an 
orbit traverses in one iteration around a 1:1 thin tube orbit). Similarly, in the vicinity of a thin orbit, ratios 
of leading frequencies will be nearly invariant, and appear as flattenings on the curve. 

The sequence of Figure 12 clearly reinforces the findings of our earlier studies of axisymmetric scale-free 
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Fig. 11. — Top: Some major thin resonant orbits, each representative of an entire family of orbits, originating 
from points in the x-z start-space inside the radius ot^hc thin tube for the logarithmic potential (L z = 0.1, 
c/a = 0.3, T = 0). Bottom: Poincare surface of section for the thick families associated with the thin 
resonant orbits shown on the left half of the figure. The islands corresponding to the higher order resonances 
occupy smaller areas. 
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Fig. 12. — Ratios of leading frequencies from the x-z start-space for a logarithmic potential at varying 
L z values L z = 0,0.1,0.2,0.3,0.4,0.5. Dashed lines represent the locations of thin tube at that angular 
momentum. Graphs for higher values of L z (up to its maximum of L z = e -1 / 2 = 0.6065) are similar to the 
graph of L z = 0.5, with the domain continually shrinking until it becomes a single point overlapping with 
the thin tube at the maximum L z . 
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potentials (Hunter et al. 1998) that resonant orbits dominate the phase space at low angular momenta, 
while the only orbits that remain at large angular momenta are tubes. It is also apparent that most of the 
resonant orbits bifurcate from the bounding outer planar orbit, while the thin tube becomes unstable only at 
low angular momenta (Hunter et al. 1998; Terzic 1998). The effectiveness of graphs of ratios of the leading 
frequencies in identifying the origin of the instabilities makes it beneficial to use them in conjunction with 
traditional Poincare surfaces of section in the analysis of phase space and stability. 

For triaxial potentials, we divide the duty of sampling the full range of orbits between two 2-D spaces, as 
outlined above, so that each picks up different types of orbits. The stationary start-space is supposed to select 
all the resonant orbits with turning points on the cquipotential surface. These will not yield any tube orbits, 
since now L z is no longer an integral of motion. The principal plane start-space is three-fold, one in each 
principal plane. It is designed to pick up tube orbits, but in the process will also pick up resonant orbits which 
pierce principal planes with velocities normal to the plane (Schwarzschild 1993). Some duplication of orbits 
is inevitable. If the only goal of the principal-plane start-space is to sample the tube orbits, restricting it to 
the area outside the thin tube orbit will indeed minimize duplication of tube orbits, without systematically 
leaving out any of them. The same is true of the stationary start-space; if its only goal is to produce the 
boxes, boxlets and other resonant orbits with zero-velocity turning points, it will indeed do so, albeit with 
some duplication, without systematically leaving out any orbits. Of course, some minor orbital families may 
be left out because of the finite resolution of the coverage of these spaces. However, at least theoretically, 
all tube orbits and orbits with zero-velocity turning points have a chance of being represented, and in the 
limit of number of points in these spaces N — ► oo, they are represented. The only remaining question then is 
whether there exist other families of orbits that this choice of start-spaces systematically leaves out, denying 
them even a theoretical chance of being represented. One cannot be absolutely definitive in answering this 
question, but all the empirical indications are that this is not the case, at least not on the level at which it 
would seriously jeopardize accurate sampling of the phase space. 
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